function diff = rep_cost_fun( m, param) 

    tmp2 = -2.*param.c.*(1-param.alpha)./param.sigma^2./(param.gamma_R_2-param.gamma_R_1).* ...
        ((1./(1-param.alpha)-param.gamma_R_1).*param.m_u.^param.gamma_R_1.*int_cost1_fun(param.m_u,param).*param.gamma_R_1./param.m_u ...
        -(1./(1-param.alpha)-param.gamma_R_2).*param.m_u.^param.gamma_R_2.*int_cost2_fun(param.m_u,param).*param.gamma_R_2./param.m_u ...    
        +(param.gamma_R_2-param.gamma_R_1).*(param.zeta+delta_fun(param.m_u,param))./param.m_u );
    tmp1=J_0_fun( param.m_h, param);

    A = (tmp1.*(param.m_u./param.m_h).^param.gamma_R_2.*param.gamma_R_2./param.m_u - tmp2) ./ ...
    ((param.m_u./param.m_h).^param.gamma_R_1.*param.gamma_R_1./param.m_u - (param.m_u./param.m_h).^param.gamma_R_2.*param.gamma_R_2./param.m_u);

    diff = - 2.*param.c.*(1-param.alpha)./param.sigma^2./(param.gamma_R_2-param.gamma_R_1).* ...
        ((1./(1-param.alpha)-param.gamma_R_1).*m.^param.gamma_R_1.*int_cost1_fun(m,param) ...
        -(1./(1-param.alpha)-param.gamma_R_2).*m.^param.gamma_R_2.*int_cost2_fun(m,param)) ... 
        + A.*(m./param.m_h).^param.gamma_R_1 ...
        - (tmp1+A).*(m./param.m_h).^param.gamma_R_2;
    diff = - diff.*(m>=param.m_h);

end